Numerical study of 3D vesicles under flow: discovery of new peculiar behaviors 
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The study of vesicles under flow, a model system for red blood cells (RBCs), is an essential step 
in understanding various intricate dynamics exhibited by RBCs in vivo and in vitro. Quantitative 
3D analyses of vesicles under flow are presented. The regions of parameters to produce tumbling 
(TB), tank-treating (TT), vacillating-breathing (VB) and even Kayaking (K) modes are determined. 
New qualitative features are found: (i) a significant widening of the VB mode region in parameter 
space upon increasing shear rate 7 and (ii) a striking robustness of period of TB and VB with 7. 
Analytical support is also provided. These findings shed new light on the dynamics of blood flow. 

PACS numbers: 87.16.D- 83.50.Ha 87.17.Jj 83.80.Lz 87.19.rh 



Introduction After nearly a century of research on 
blood, understanding of the basic blood flow mecha- 
nisms at the cellular level (at the scale of red blood cells, 
platelets, etc.) is still an open issue. Blood is a complex 
fluid and the description of its flow properties escapes 
the traditional Navier-Stokes law known for simple fluids 
(e.g. water). 

To date blood flow has been described by means of 
phenomenological continuum models that require many 
assumptions which are difficult both to justify and to 
validate Furthermore, the modern view of complex 
fluids has highlighted the intimate coupling between the 
dynamics of the microscopic suspended entities (RBCs 
in the example of blood) and the flow at the global scale 
[2] • This implies that the macroscopic law of blood flow 
should carry information on the microscopic scale, such 
as the orientation of the cells, their individual and col- 
lective dynamics, their local concentration (hematocrit), 
and so on. Apart from a dilute suspension where a rhe- 
ological law can be extracted analytically [3], a complete 
understanding of blood flow should ultimately emerge 
from a numerical study. 

Computational approaches are, however, challenging 
due to the free boundary character of the RBCs (the 
shape is not known a priori) which is fixed via a subtle 
interplay between the local flow and the different internal 
modes of the RBCs (membrane bending, shear elasticity). 

Under simple shear flow, RBCs exhibit tank-treading 
(TT), tumbling (TB), vacillating-breathing (VB) (aka 
swinging or trembling) and kayaking (K) (see movies 
of the four modes [3]). These dynamics should have an 
impact on the macroscopic flow. These modes are shared 
both by vesicles @, H, [fl 0, H ED, GS E Q EE G3, E3| 
(made of a pure bilayer phospholipid membrane) and 
RBCs flg| . which, apart from having the same type of 
phospholipid membrane, are endowed with a cytoskele- 
ton (a cross-linked network of proteins lying underneath 
the RBC membrane). A systematic analysis of these dy- 



Figure 1: A typical vesicle obtained by simulations. The white 
line indicates the main axis direction. The color code on the 
vesicle corresponds to the value of the local curvature, which 
is large at the tips and low in the central region. 



namics constitutes an essential step that is explored in 
this Letter. 

We consider a model system, which is a phospholipid 
vesicle (FigureQ}. We thus ignore in this first exploration 
the effect of the cytoskeleton. This reduction is essential 
to identify the basic underlying mechanisms. . 

The major results of this study are: (i) We find that 
the width of the region of parameters where the VB mode 
is manifested significantly widens upon increasing shear 
rate. This behavior was not captured by any of the pre- 
vious theories or simulations, (ii) We find that the period 
of oscillation in the TT and TB regime is practically in- 
dependent of the shear rate 7, in marked contrast with 
a previous numerical study using multiparticle collision 
dynamics [l6| • (hi) We develop an analytical theory that 
accounts for these two facts, (iv) We demonstrate that 
the Keller-Skalak theory (which assumes a fixed shape 
for the RBCs, and is often used as a basis in experiments 
and theoretical studies) is an excellent framework for TT 
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and TB description when shape deformability is negli- 
gible (small 7 or large A, the ratio of internal over the 
external viscosity). However, major deviations are re- 
vealed when 7 is large enough and/or A is not too large, 
(v) We report on a transition between tumbling and the 
kayaking phase, during which the main axis of the vesicle 
rotates around an axis perpendicular to the shear plane. 
Strikingly, kayaking seems to occur even if the vesicle is 
initially forced to be in the shear plane. 

The Model: The Stokes equations (experiments have 
so far explored the limit of small Reynolds numbers) can 
be formally solved using the Boundary Integral (BI) for- 
malism (see for example (2l|) which yields 

1m v mem( r ) = Vout^ shear + J mem G( r ~ r )^mem( r ) dr 

+ (Vin - Vout) J mem v mem (r ).K(t - r').n(r') dr 

where v mem (r) is the local velocity field at the mem- 
brane, v s hear = 7£/x is the externally applied Couette 
flow (with 7 the shear rate), r\ m = (??,„ + Vout)/ '2 (rj in 
and rjou stand for viscosities of the internal and external 
fluid, respectively), J mem is an integration over the mem- 
brane and G(r — r ) and K(r — r ) are the Green tensors 
defined as: 
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Results: The phase diagram of various dynamics (TT, 
TB and VB) of a vesicle and RBC in a Couette flow re- 
cently motivated a large number of theoretical and ex- 
perimental studies @, EH E3, G3 HO]. This is the first 
focus of our study. Let us introduce appropriate dimen- 
sionless numbers. Shape deformability can be measured 
by the dimensionless number 



A vesicle that is subjected to a flow undergoes a shape 
transformation that is limited by bending modes. The 
reaction bending force of a vesicle on the fluid is given 
by the Helfrich force Q 

fcun>(r) = -« \c{y) {c(r) 2 - 4 5 (r)} + A 2D c{r) n(r) 

(2) 

where K is the bending modulus, c(r) is the local mean 
curvature (c(r) = ci(r) + czir), c\ and c 2 are the two 
principal curvatures at point r of the membrane), g(r) 
is the Gaussian curvature (<? = C1C2), A 2 jj is the (sur- 
face) Laplace-Beltrami operator on the membrane, and 
n is the outward unit normal vector. The other contribu- 
tion to the membrane force follows from local membrane 
incompressibility: 

f tms (r) = T[C(r)c(r)n(r) + V 2D ((v)} 

where V 2 d is the surface gradient operator and £(r) is 
a local dimensionless Lagrange multiplier that enforces 
membrane inextensibility (the constant T has a dimen- 
sion of energy per unit surface). Local membrane inex- 
tensibility sets severe limitations on the numerics, both 
on the time step and the precision of the results. Care- 
lessness (for example, imposing a small enough toler- 
ance of a few % on area) may induce spurious results. 
The scheme used in this study ensures an area variation 
lower than 5.10 _4 %. The precise numerical analysis is 
a formidable task and technical details will be reported 
elsewhere. Instead, we focus here on the major outcomes. 



where R is the typical size of the vesicle. Vesicles exhibit 
various equilibrium shapes (i.e., in the absence of flow) 
[13], depending on their reduced volume r = S 3 ^ 2 , 
where S is the area of the membrane and V the inter- 
nal volume, t is the second dimensionless parameter. 
For t > 0.652, the shape is prolate (one long revolu- 
tion axis) In the range 0.591 < r < 0.652, the 
shape is oblate (one small revolution axis; the shape is 
biconcave, known also for RBCs). We found that, for 
C K > 1 (which is quite easily reachable experimentally), 
the oblate branch is suppressed in a Couette flow so that 
only prolate shapes prevail (see Fig[T]). The last dimen- 
sionless parameter is A = rji n /r] ou t. While our study 
can be performed at different r, we will consider a re- 
duced volume close enough to unity. This will allow us 
to compare the theoretical approaches considering the 
quasi-spherical limit l3, 14, 15, 0, [Tt]]. We have cho- 
sen the value r = 0.95 and have explored the effects of 
the two other parameters C K and A. The resulting phase 
diagram is represented in FigfSl 

A first important feature found here is that the bound- 
aries of the phase diagram are strongly underestimated 
in previous analytical theories 15J, ll7|]. For example, for 
C K ~ 0.1, the bifurcation from TT to TB occurs here for 
A ~ 8, while analytical theories predicted A ~ 4 (Figf2j. 
This is astonishing because deviation from a sphere is 
only about 5% and a perturbative scheme is expected to 
make sense. Thus, we have attempted to understand this 
peculiar behavior by revising previous analytical theories. 
The key ingredient is that the next order terms in an ex- 
pansion in powers of excess area (relative to a sphere) are 
decisive, however small the deviation from a sphere. We 
have identified the fact that this is the result of a singular 
behavior of the expansion scheme. Details of the analyt- 
ical theory will be given elsewhere [24|. The outcome of 
this analytical theory is presented in Figure [21 revealing 
excellent agreement with the full numerical simulation. 

Previous analytical [3, EH] as well as numerical calcu- 
lations [H] (based on dissipative particle dynamics) and 
experiments [iji] reported that the band of existence of 
the VB mode saturates with C K above a value of the or- 
der of C K ~ 0.5. This is in good agreement with the 
present numerical study as long as 0.5 < C K < 2. Be- 
yond a value of the order of C K = 2, the VB band exhibits 
a sudden ample widening as shown in Fig. [2j This ef- 
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Figure 2: Phase diagram for r = 0.95. The symbols corre- 
spond to numerical results and the solid line to the analytical. 
The Flow alignment states correspond to a VB mode with a 
very small amplitude (see Fig[5j) . These states are thus inter- 
mediate between VB and TT. The K phase is visible at the 
upper right corner of the diagram. 



feet highlights the nontrivial character of dynamics due 
to shape deformation. A careful theoretical analysis of 
this phenomenon has led us to the discovery that not 
only is the fourth order harmonic strongly coupled to 
the second one, but it also acquires strong activity on 
increasing C K . This mode, which is damped at low de- 
formability, becomes active (i.e. it is excited) at larger 
deformability (larger C K ). As a consequence, the VB 
regime is promoted further and further, leading to the 
VB band widening. We have found that this peculiar 
behavior is captured by the new theory that implements 
fourth order harmonics 

Let us now analyze the behavior of the tumbling angle 
9(t). It is convenient to represent a phase portrait in the 
plane (8,9). This will also allow us to shed light on the 
limit of applicability of the KS theory, which is often used 
as a basis in experiments [H|,[I3] and in numerical simula- 
tions [H]. We find that, for C K -> (we choose C K = 0.1 
as an example), the simple relation 9/j = A + B cos(29) 
predicted by the KS theory is in excellent agreement with 
the BI simulations. 

At larger deformabilities, the situation is more com- 
plex and requires a careful analysis. We first focus on the 
value A = 20. The result is reported in Figj3] (shifted by 
+0.5 upward for clarity). All data for 1 < C K < 10 are 
nearly superimposed, reflecting a good agreement with 
the KS theory even though, for C K > 5, the vesicle is 
actually kayaking [IB]. 0(t) represents in this case the 
tumbling of the projection of the main axis in the shear 
plane. In contrast, at smaller values of A, (A = 16 in 
FigEJ) , major deviations from the KS theory are man- 
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Figure 3: 6 as a function of 9 for A = 16, A = 20 (shifted 
upwards by +0.5), for various C K . The dashed line cor- 
responds to the KS theory. The full black line is a fit of 
the data points forA = 16 and C K = 8 with the func- 
tion A + Bcos(260 + Ccos(40) (A = 0.543, B = 0.462 and 
C — 0.049), showing that higher harmonics play an impor- 
tant role at large deformabilities. 



ifested. In particular, we observe the excitation of the 
fourth order harmonic, represented by cos(4£>) in the fig- 
ure when Cfc is increased up to 8. Note that, for C K = 10 
, other higher order harmonics are excited as well, which 
is a precursor to the TB-VB or the TB-K bifurcation. 

We have analyzed several other physical quantities, 
such as evolution and possible exchanges of the two main 
axes lying in the shear plane, but we shall focus only on 
the behavior of the period (scaled by the shear rate) in 
this brief exploration. The results are shown in FigHl 
A remarkable fact is the quasi-robust character of the 
rescaled period as a function of the deformation regard- 
less of the dynamical mode (TB, VB or K) (only minor 
variations by few % are manifested over a range of a 
decade in shear rate). This finding is not intuitive in as 
much as deformation (at large enough C K ) plays an es- 
sential role (strong deviation from the KS theory). This 
result runs contrary to a previous numerical analysis [l6| 
based on dissipative particle dynamics. The ratio of their 
period variation and ours attains a factor 100. We do not 
have an explanation about the origin of this major dis- 
crepancy. We have checked, in light of the new analytical 
theory, that the same quasi-independence on C K is also 
manifested. 

Finally, it is worth emphasizing that the transient dy- 
namical regime can prove to be very slow, especially at 
large values of C K . It is necessary to ascertain the es- 
tablishment of a given regime on the scale of the longest 
time. For large enough A, the appropriate (slow) time 
scale is r K = rji n R? / n; motion is limited by the more 
viscous internal fluid. For C K = 10, this may typically 
represent a hundred cycles in the TB regime. Thus, it is 
necessary to make use of appropriate and stable numeri- 
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Figure 5: Long time relaxation of the vesicle close to the 
TT-VB line. The insert shows a log plot of the vacillation 
amplitude. 



cal schemes. Details will be published elsewhere. Here, it 
suffices to provide an illustration exhibiting the slow re- 
laxation (FigEJ. The VB amplitude close to the TT-VB 
line (A = 9, C K — 8) is shown. We find an exponential 
relaxation of the VB amplitude (the insert of FigE] pro- 
vides a log plot of the amplitude). The vesicle reaches 
its steady state only after a characteristic time of about 
30r]i n R 3 1 k, corresponding to 125 oscillations. The same 
feature is observed for the TB-K transition. 
Conclusion: We have reported on a quantitative numer- 
ical analysis in three dimensions that led us to identify 
new features of vesicle dynamics, which were not revealed 
in prior analytical, numerical, or experimental studies. 



A theory has allowed us to unearth a subtle role played 
by higher order terms. The theoretical results agree re- 
markably well with the numerical ones. This work is a 
first essential step toward a systematic study of RBC dy- 
namics and blood rheology. An obvious limitation of our 
work is the absence of shear elasticity associated with 
the cytoskeleton. An analysis that includes this factor 
constitutes the next natural step. 

CM. and A.F. Acknowledge financial support from 
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